#ifndef PERLIN_H
#define PERLIN_H

#include <cmath>

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif



double Cubic_Interpolate(double v0, double v1, double v2, double v3, double x)
{
	double P = (v3 - v2) - (v0 - v1);
	double Q = (v0 - v1) - P;
	double R = v2 - v0;
	double S = v1;

	return P*pow(x,3) + Q*pow(x,2) + R*x + S;
}

double Cosine_Interpolate(double a, double b, double x)
{
  double ft = x * 3.1415927;
	double f = (1 - cos(ft)) * .5;

	return  a*(1-f) + b*f;
}

double Noise(int x, int y)
{
  int n = x+y*57;
  n = (n<<13)^n;
  return (1.0 - ( (n*(n*n*17731+789221) + 1376312589) & 0x7fffffff) / 1073741824.0);
}

double SmoothedNoise1(double x, double y)
{
  double corners = ( Noise(x-1, y-1)+Noise(x+1, y-1)+Noise(x-1, y+1)+Noise(x+1, y+1) ) / 16;
  double sides   = ( Noise(x-1, y)  +Noise(x+1, y)  +Noise(x, y-1)  +Noise(x, y+1) ) /  8;
  double center  =  Noise(x, y) / 4;
  return corners + sides + center;
}

double InterpolatedNoise_1(double x, double y)
{
  int integer_X    = int(x);
  double fractional_X = x - integer_X;

  int integer_Y    = int(y);
  double fractional_Y = y - integer_Y;

  double v1 = SmoothedNoise1(integer_X,     integer_Y);
  double v2 = SmoothedNoise1(integer_X + 1, integer_Y);
  double v3 = SmoothedNoise1(integer_X,     integer_Y + 1);
  double v4 = SmoothedNoise1(integer_X + 1, integer_Y + 1);

  double i1 = Cosine_Interpolate(v1 , v2 , fractional_X);
  double i2 = Cosine_Interpolate(v3 , v4 , fractional_X);

  return Cosine_Interpolate(i1 , i2 , fractional_Y);
}

double PerlinNoise_2D(double x, double y, double per, double Oct)
{
  double total = 0;
  double p = per;

  for(int i=0; i<Oct; i++)
  {
    double frequency = pow(2.0,i);
    double amplitude = pow(p,i);

    total = total + InterpolatedNoise_1(x * frequency, y * frequency) * amplitude;
  }

  return total;
}

double *baseNoise; //an array for our uniform noise
double *smooth;
double *perlinNoise;
int width, height;

void generateUniformNoise()
{
	for(int i = 0; i < width; i++)
		for(int j = 0; j < height; j++)
			//random float between 0 and 1
			baseNoise[j*width + i] = unif();
}

double interpolate(double x0, double x1, double alpha) //alpha lies between 0 and 1
{
	return (1.0 - alpha) * x0 + alpha * x1;
}

void generateSmoothNoise(int k)
{
	int samplePeriod = 1 << k; // calculates 2 ^ k
	double sampleFrequency = 1.0f / samplePeriod;
	for(int i = 0; i < width; i++)
	{
		//calculate the horizontal sampling indices
		int sample_i0 = (i / samplePeriod) * samplePeriod;
		int sample_i1 = (sample_i0 + samplePeriod) % width; //wrap around
		double horizontal_blend = (i - sample_i0) * sampleFrequency;
		for(int j = 0; j < height; j++)
		{
			//calculate the vertical sampling indices
			int sample_j0 = (j / samplePeriod) * samplePeriod;
			int sample_j1 = (sample_j0 + samplePeriod) % height; //wrap around
			double vertical_blend = (j - sample_j0) * sampleFrequency;
			//blend the top two corners
			double top = interpolate(baseNoise[sample_i0 +width*sample_j0],
                               baseNoise[sample_i1 +width*sample_j0],
                               horizontal_blend);
			//blend the bottom two corners
			double bottom = interpolate(baseNoise[sample_i0 +width*sample_j1],
                                  baseNoise[sample_i1 + width*sample_j1],
                                  horizontal_blend);
			//final blend
			smooth[j*width + i] = interpolate(top, bottom, vertical_blend);
		}
	}
}

void generatePerlin(int octaveCount)
{
	double persistance = 0.5f;
	//generate smooth noise

	double amplitude = 1.0f;
	double totalAmplitude = 0.0f;
  double weight = 0.5;
	//blend noise together
	for(int k = 0; k < octaveCount; k++)
	{
    generateSmoothNoise(k);
		amplitude *= persistance;
		totalAmplitude += amplitude;

		for(int i = 0; i < width; i++)
			for(int j = 0; j < height; j++)
				perlinNoise[i +width*j] += smooth[i +width*j] * weight;
    weight *= 0.5;
	}
	//normalisation
	for(int k = 0; k < octaveCount; k++)
	{
    for(int i = 0; i < width; i++)
			for(int j = 0; j < height; j++)
        perlinNoise[i +width*j] /= totalAmplitude;
	}
}

void setup(int w, int h)
{
  width = w;
  height = h;
  baseNoise = (double *)malloc(sizeof(double)*width*height);
  smooth = (double *)malloc(sizeof(double)*width*height);
  perlinNoise = (double *)malloc(sizeof(double)*width*height);
  for(int i = 0; i < width; i++)
    for(int j = 0; j < height; j++)
      perlinNoise[i +width*j] = 0;
  generateUniformNoise();
  generatePerlin(4);
}

#endif // PERLIN_H
